Modeling near-field radiative heat transfer from sharp objects using a general 3d 

numerical scattering technique 
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We examine the non-equilibrium radiative heat transfer between a plate and finite cylinders and 
cones, making the first accurate theoretical predictions for the total heat transfer and the spatial heat 
flux profile for three-dimensional compact objects including corners or tips. We find qualitatively 
different scaling laws for conical shapes at small separations, and in contrast to a fiat/slightly- 
curved object, a sharp cone exhibits a local minimum in the spatially resolved heat fiux directly 
below the tip. The method we develop, in which a scattering-theory formulation of thermal transfer 
is combined with a boundary-element method for computing scattering matrices, can be applied to 
three-dimensional objects of arbitrary shape. 



Introduction: We make the first accurate theoretical 
predictions for near-field thermal transfer from 3d com- 
pact objects of arbitrary shape (including corners or tips) 
to a dielectric substrate. Our work is motivated by stud- 
ies of non-contact thermal writing with a hot, sharp 
object [U [2]- Theory has predicted [31 E] and experi- 
ments have confirmed [51 [6] that radiative heat transfer 
between two bodies at different temperatures is greatly 
enhanced as their separation is reduced to sub-micron 
scales, due to contributions from evanescent waves. Un- 
til the last few years, the only rigorous theoretical re- 
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FIG. I: Total thermal transfer between a silica plate and 
doped silicon objects of various shapes. The plate is semi- 
infinite, and the objects all have height equal to 1 /xm in the z- 
direction. The plate/environment temperature is Tp = 300 K 
and the objects are at temperature Ta = 600 K. Red dots 
denote results with the sphere scattering matrix determined 
analytically, the only case in which we have an analytic solu- 
tion for the scattering matrix. 



suits for thermal transfer concerned parallel plates; how- 
ever, very recently rigorous theoretical predictions for 
sphere-sphere [7] and sphere-plate [HI [5] geometries as 
well as general formalisms for planar structures [ID] and 
arbitrary shapes [H |TT] have been presented. Neverthe- 
less, such techniques were previously implemented only 
when analytic expressions for the scattering matrices 
were known (e.g., spheres and plates in 3d). As an al- 
ternative, stochastic finite-difference time-domain meth- 
ods have been used to examine heat transfer for periodic 
structures [V2\ . but this method is not computationally 
well-suited for compact objects in three dimensions. Our 
technique extends the formalism of Ref. 151 directly to ar- 
bitrary compact objects. To do this, we use a boundary- 
element method in which the object is described by a 
generic surface mesh jl3| . We then numerically compute 
the scattering matrices of this object in a multipole ba- 
sis; for our study, we employ a cylindrical-wave basis. 
Unlike the usual spherical-wave basis, this allows us to 
concentrate our resolution on the surfaces adjacent to the 
substrate, but requires a new quadrature approach to dis- 
cretize the scattering matrix. In addition to sphere-plate 
heat transfer, we study both cylinder-plate and cone- 
plate configurations (see sketch in Fig. fTl), for which no 
known analytic solution exists. Our results exhibit clear 
scaling laws for the total heat transfer that distinguish 
locally flat structures (e.g., cylinders and spheres), from 
locally sharp structures (cones). In addition, we study 
the spatial distribution of heat-flux over the substrate, 
a topic that has been treated previously using a point- 
dipole approximation for the heat source [14] . Our results 
show that the heat flux pattern depends strongly on the 
shape of the tip. Cones in particular have a flux pattern 
exhibiting an unusual feature: a local minimum in the 
heat flux directly below the tip, which we can explain 
with a modified dipole picture. 

Method: In our setup, an object A at (local) tempera- 
ture Ta faces a dielectric plate P at temperature Tp, in 
an environment E that is also at temperature Tp. We use 



the framework of Rytov's theory [T5] , in which all sources 
emit radiation independently. The full non-equilibrium 
Poynting flux can be computed with radiative sources 
from P and E only, as the flux from A at temperature 
Ta must equal the flux from P and E at temperature 
Ta (with opposite sign), due to detailed balance p^ . 
To compute the power flux, we first compute the non- 
equilibrium electric field correlator (E(x) (g)E*(x'))j due 
to radiation from j = P,E for general x 7^ x' (the 
Poynting fiux will be obtained at the end by taking 
lima;/_j.a; Va:' (E ig) E'*)) which is expressed as an integral 
of the general form: 
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where 9 = cj* (exp (/iw/fesT) - 1)"^ [H US], h is 
Planck's constant and k^ the Boltzmann constant. Un- 
less otherwise noted, we consider each frequency u sepa- 
rately and drop the uj subscript below. 

The correlator takes on a simple form in an orthogo- 
nal basis Eq((u;; x) for the field degrees of freedom (in our 
case, these will be cylindrical waves in the ±z direction), 
indexed by a (discrete or continuous) index a, and rep- 
resent the correlator as a matrix D. In matrix notation 
(with implied summation over repeated indices): 



(E(x)®E*(x')) 



_„E„,(x)®E^(x'). (1) 



D = Dp +De due to statistical independence of the ther- 
mal fluctuations, where Dp/^ involve sources only from P 
I E. The correlators Dp/£; are obtained from the "unper- 
turbed" correlators Dp;^;; Dp involves the plate sources 
without A and D^ involves the environment sources with 
neither A nor P present. The D^ are known analyti- 
cally (see below), and the full correlators Dj can be de- 
termined from them by use of the Lippmann-Schwinger 
equation [HI [IT'. In our notation: 
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The Oj are matrices that describe the scattering of in- 
coming and outgoing fields with the allowance for sources 
in between the objects, described explicitly in [TH]. 
These are constructed from the more conventional incom- 
ing/outgoing scattering matrices Fp/^ [l71[T9] for objects 
P and A individually. As object P is a plate, Fp is known 
analytically. However, ^ a cannot be determined analyt- 
ically for a general object A. Instead, the computation 
of the scattering matrix elements is accomplished via a 
boundary-element method [13], described below. The z- 
component of the Poynting fiux at position x, S'x, and the 
total power fiux St through the z = plane can both be 
expressed as operator traces: ^x/t — Re Tr [Sx/tDt] , 
with (Sx)„,,„ = -^z • [E„(x)x (VxE„,(x))*] and 
^T)a',a gi'^en below. 

We employ a cylindrical-wave basis of fields 
Es,m,fc ,p(x) in which the waves (also known as 



Bessel beams) propagate in the ±z direction [50!. The 
variable s = ± refers to the direction of propagation; m 
is the (integer) angular moment of the field, < fcp < 00 
the radial wavevector, and p = Af , N the polarization. 
The composite index in this case is a = {s,m,kp^j)\. 
This basis is especially well-suited to the case considered 
here in which objects have rotational symmetry about 
the z-axis, as different values of m are decoupled. 

To compute the elements of F^, we use a boundary- 
element method (BEM) [T31 [H]. In this framework, 
the surface of object A is discretized into a mesh; 
our numerical method then computes the induced cur- 
rents from an incident multipole field Eq,(x) (here a = 
{s, m, fcp,^}). The multipole moments of this current dis- 
tribution are then computed in a straightforward man- 
ner [5D], which yield the scattering matrix F^ W^.- Be- 
cause the cylindrical-wave basis distinguishes between 
waves in the ±z direction (unlike a spherical wave ba- 
sis) , and because the near-field thermal transport mostly 
depends on reflections from adjacent surfaces, we are able 
to concentrate most of our BEM mesh resolution on the 
part of the surface of A nearby the plate, greatly improv- 
ing computational efficiency. For example, in the mesh 
for a cone below we use ~ 250 times more resolution at 
the tip than at the base. 

One complication of cylindrical multipoles is that fcp is 
a continuous index and matrix multiplication is turned 
to integration. For computational purposes, this inte- 
gration must be approximated as a discrete sum by nu- 
merical quadrature. We approximate the integral over 
fcp using a Gaussian quadrature scheme [55] for high ac- 
curacy. For example, consider the scattering matrix V a 
of object A\ its action on an incident electric field can 
be discretized as (for simplicity, summation over m and 

p is suppressed): F^E^ . 
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Z]j=o'^j C^A)j.i^fep,j where the sets {w^, fcp^} form a 
set of one-dimensional quadrature weights and points, 
respectively, and (F^) ^ = (Fa)j. , ^ , are the elements 
of the continuous scattering matrix. 

The analytic expression for the non-equilibrium elec- 
tric field correlator of a plate at temperature Tp and 
environment at T = expressed in the planewave basis 
is well-known [T^ [TB] . Since there is a standard identity 
relating planewaves to cylindrical waves, it is a simple 
exercise to re-express this correlator in the basis of cylin- 
drical multipoles [50] : 
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Here r^ p are the Fresnel coefficients for a dielectric 
plate, Xp(e) = 1 for kp < uj {kp > w) and zero other- 
wise, q — \/oj^ ~ fcp, and 6ij is the Kronecker (Dirac) 
delta function on discrete (continuous) indices; the Ss.+ 



reflects the fact that only waves propagating in the +z 
direction are emitted by the plate. The expression for 
the environment correlator D^ is given by the same ex- 
pression as Dp with r = and 5s^+ replaced with 5s ~- 
Finally, the matrix elements for the total power flux are 

For the surface meshes, we use approximately 2,500 
panels (discretized surface elements) to get 1% conver- 
gence, with the panels highly concentrated on the area 
of the objects nearest to the plate. We retain angular mo- 
ments up to \m\ = 10, and for each m we perform the w 
and kp integrations using 28 and 48 Gaussian quadrature 
points, respectively. For our study, object A is composed 
of doped silicon while the substrate B is silica. For the 
doped silicon dispersion we use a standard Drude-Lorentz 
model 1231 with a dopant density of 1.4 x lO^^cm"'^, while 
for silica we use measured optical data [5]. 

Results: Figure [T] shows the geometry-dependence of 
the total heat transfer rate between different compact ob- 
jects and a dielectric plate, over surface-surface separa- 
tions z from several microns down to 20 nm. In addition 
to the expected near-field enhancement, we observe sev- 
eral crossings as, e.g., the broader surface area of the 
R — 0.5 //m radius sphere competes with the smaller 
but flatter surface of the d = 0.4 /xm diameter cylin- 
der. For smaller z, the ratio of the transfer between the 
d — 0.4 /.tm and d — 0.2 /im cylinders approaches the 
ratio of their surface areas (within 6% at z = 20 nm) , 
as would be expected from a proximity approximation 
(PA) [S1I21]- The sphere-plate exhibits the 1/z power law 
as predicted by PA IB |1] to within 10% for z < 0.1 /im, 
while the cylinder-plate exhibits agreement to within ap- 
proximately 10% over this range using a PA based on 
the integral of the plate-plate heat transfer rate over the 
cylinder front face and vertical sidewalls. The contri- 
bution from the sidewalls can be ignored (leading to a 
~ 1/z^ transfer rate [J) for z/d < 0.01. In contrast to 
the sphere and cylinders, the cones do not seem to be 
asymptoting to a power law, and may even have a loga- 
rithmic divergence as z — > 0, a fact which we attribute to 
the scale-invariance of the plate-cone configuration when 
z ^ 1/xm and z <C hc/ksT (the latter eliminating ma- 
terial dispersion effects). To check the accuracy of our 
numerical scattering method, we also plot the results for 
the sphere where F^ is calculated semi-analytically [S], 
shown as red dots, which agrees to within 1%. 

For thermal writing applications, an important factor 
to consider is not only the total power delivered to the 
plate, but also the spatial extent over which this delivery 
occurs. In order to examine this, we envision a scenario in 
which a critical magnitude of the z-directed Poynting flux 
is required in order for some change to occur on the plate, 
for example, the patterning of a thermal mask for later 
etching [5] . Figure [2] plots the Poynting flux at x = as 
a function of z, which will tell us how far away the object 
must be before it can effect this patterning. The cylin- 
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FIG. 2: Poynting flux at the origin for the geometries of Fig. [T] 
with plate/environment temperature Tp = 300 K and ob- 
ject temperature Ta = 600 K, using the single-polarization 
approximation (SPA). The magenta line denotes the sphere- 
plate without the SPA, and the horizontal dashed line denotes 
the threshold used for the crossectional flux profiles of Fig. [3] 



ders and spheres converge to the same ~ 1/z^ profile for 
small z (as expected from a PA), whereas the cones all 
follow 1/z^ profiles with different coefficients. This 1/z^ 
dependence follows from the scale-invariance of the scat- 
tering problem for small z, combined with the fact that 
there is a 1/z cutoff in the range of kp that contributes 
to the transfer, so that the total number of modes that 

1 It 

contribute is proportional to /^ dkpkp ^ 1/z^. In this 
calculation we have found that the result is dominated 
by the N polarization (E _L z), mirroring similar phe- 
nomena in other near-field cases [5S], that results from 
the behavior of the Fresnel coefficients for high kp. This 
is fortunate because we have found that the M contri- 
bution to the Poynting ffux requires much higher mesh 
resolution to converge. To check this single-polarization 
approximation (SPA) for a sphere we also plot the full 
results, finding that the error from the SPA is < 20% at 
the largest z, decaying to < 10% at smaller z; SPA for a 
cone is discussed below. 

Figure |3] plots the Poynting flux as a function of 
X showing the heat transfer profile. For each object, 
we chose z to have the same x = Q Poynting flux of 
10~^(27r)^/ic//im^ (horizontal dashed line in Fig. l2|, cor- 
responding to a sphere-plate separation of w 200 nm. The 
cylinders and 120° cone all reach this threshold at com- 
parable separations, whereas the 90° cone is at less than 
half the separation, and the 40° cone does not even reach 
this threshold within the range considered. 

Fixing the peak Poynting flux to 10^'^(27r)^ft,c//im'', 
in Fig. |3] we plot the Poynting flux profiles for these 
shapes as a function of x. The widths for the cylinders are 
narrower than the sphere, implying that the cylinders can 
write higher spatial resolution. Surprisingly, the cones do 
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FIG. 3: Spatially-resolved heat flux profiles at the substrate 
surface. z is chosen to fix Poynting flux at a; = at 
10"^(27r)^/ic/^m^. 
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FIG. 4: Spatially-resolved heat flux profiles (arbitrary units) 
at the substrate surface for three cones at a single frequency 
cj = 0.3066(27rc/^m) and fixed z = 70 nm. The profiles 
are normalized so that their maximal value is equal to 1. 
These profiles are computed without the SPA using much 
finer meshes. For comparison, the profile for a cylinder of 
radius d = 400 nm (using the SPA) is shown as well. 



not exhibit this simple behavior. Rather, the Poynting 
flux profiles for the two cones are noK-monotonic in x, 
with a local minimum at a; = 0. The degree of non- 
monotonicity appears to increase as the cone becomes 
sharper. 

Before attempting to explain this effect, we must first 
recall that the results of Fig. [2] and Fig. |3] relied on the 
SPA; although we know this approximation to work well 
for flat or smoothly curved bodies, it is not obvious that 
it applies equally well to the cone. To confirm this result 
without this approximation, we must go to a much denser 
mesh near the cone tip to ensure mesh convergence; for 
this, we form a mesh using approximately 12,000 panels 
for these cones. We have observed that for the shapes and 
separations of interest here, the Poynting flux profiles at 



all relevant frequencies have very similar shape, and are 
simply scaled by a frequency- and material-dependent 
weight. Therefore, it is sufficient to consider a single 
frequency, which we pick to be w = G.3066(27rc//im). 
The resulting Poynting flux profiles for all three cones 
at a fixed z — 0.1 /im are shown in Fig. [4J for ease of 
comparison, all curves are scaled to have a maximum of 
1. We also show the d = 0.4 //m cylinder (using the SPA) 
for comparison. The dip at x = is less pronounced for 
the exact curves than for the SPA; in fact, the dip has 
vanished for 9 = 120°. However, it is still present for 
9 = 90° and is very prominent for 9 — 40°, where the 
Poynting flux at a; = is less than half of its peak value. 
Therefore, we conclude that this effect is not a result of 
our approximations. 

We believe the explanation for the dip in the Poynting 
flux is that as the cone tip becomes sharper, its radiation 
pattern approaches that of a dipole with axis normal to 
the plate, which has zero Poynting flux at a; = 0. This 
explanation predicts that a very thin cylinder with d <^ z 
should also have a dip in the Poynting flux at x = 0, 
which we have also confirmed numerically. 
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